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We  consider  the  effects  of  noise  on  a  model  of  epidemic  outbreaks,  where  the  outbreaks  appear  randomly. 
Using  a  constructive  transition  approach  that  predicts  large  outbreaks  prior  to  their  occurrence,  we  derive  an 
adaptive  control  scheme  that  prevents  large  outbreaks  from  occurring.  The  theory  is  applicable  to  a  wide  range 
of  stochastic  processes  with  underlying  deterministic  structure. 
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I.  INTRODUCTION 

Recently,  there  has  been  much  research  of  steady  state 
epidemics  in  random  populations  [1]  and  their  control  [2]. 
Nonequilibrium  diseases,  in  contrast,  are  those  diseases  ex¬ 
hibiting  outbreaks  that  fluctuate  in  time.  Childhood  [3,4]  and 
tropical  diseases  [5,6]  are  a  few  examples  of  outbreaks  hav¬ 
ing  strong  aimual  oscillations  with  random  amplitude.  In 
modeling  the  annual  incidence  of  infections,  random  compo¬ 
nents  from  the  environment  and/or  populations  play  a  sig¬ 
nificant  role  [7,8].  While  excellent  data  from  seasonally  fluc¬ 
tuating  diseases  illustrate  strong  annual  oscillations  with 
random  peak  outbreaks  in  the  infections  [4,9],  models  and 
data  analysis  reveal  that  outbreaks  stem  fi-om  stochastic  per¬ 
turbations  in  either  population  or  epidemic  parameters,  mak¬ 
ing  deterministic  prediction  difficult. 

Predictability  of  seasonally  driven  diseases  that  are  sto¬ 
chastic  is  necessary  for  the  application  of  methods  to  sup¬ 
press  future  outbreaks.  Many  vaccine  schemes  are  available 
for  equilibrium  diseases  [3,10],  but  in  the  case  of  nonequi¬ 
librium  outbreaks,  current  methods  may  enhance  outbreaks 
or  fail  to  achieve  their  goals  [11,12].  (Similar  problems  arise 
in  the  large  fluctuation  theory  of  stochastic  dynamical  sys¬ 
tems  [13].)  Other  methods  pulse  the  population  without  sam¬ 
pling  for  prediction  [14],  or  they  rely  on  reducing  spread  via 
mean  threshold  reduction  [3].  To  address  the  problem  of  sup¬ 
pressing  outbreaks  in  stochastic  epidemics,  we  apply  a  math¬ 
ematical  method  [15]  to  a  stochastic  model  to  predict  out¬ 
breaks  before  they  occur,  and  then  adapt  a  vaccine  strategy 
which  prevents  the  outbreak  firom  occurring.  The  theory  ex¬ 
ploits  a  transition  probability  description  from  small  ampli¬ 
tude  incidence  to  outbreak  dynamics,  and  generates  a  region 
of  high  probability  transport  of  the  most  sensitive  regions  to 
stochastic  effects.  Moreover,  it  allows  us  to  monitor  regions 
of  stochastic  dynamics  that  have  a  high  probability  of  pre¬ 
ceding  a  large  outbreak,  which  in  turn  leads  to  a  design  of  a 
vaccine  control  strategy  to  suppress  outbreaks.  We  thus  argue 
a  general  simple,  but  effective,  control  technique  that  takes 
advantage  of  complicated  interactions  of  determinism  and 


noise.  The  techniques  introduced  here  may  also  be  applied  to 
general  stochastic  nonautonomous  systems  of  the  form 

^  =  G{x,t)+  7j{t),  (1) 

at 

where  G{x,t)=G{x,t+\),  and  the  noise  is  added  periodically 
with  the  period  of  drive,  i.e., 

ri{t)=  rf„lS{t-n),n  =  \,2,  ,  (2) 

A  is  the  Dirac  delta  function,  and  r)„  is  now  a  discrete  ran¬ 
dom  variable.  The  form  of  Eq.  (1)  allows  us  to  consider  the 
dynamics  as  a  discrete-time  constantly  perturbed  stochastic 
dynamical  system. 

n.  A  STOCHASTIC  EPIDEMIC  MODEL 

A  standard  system  used  to  study  and  predict  the  stochastic 
dynamics  of  disease  epidemics  is  based  on  a  simplified  re¬ 
duced  version  of  the  well-known  SEIR  (defined  below)  com- 
partmental  model  [7,9,16],  known  as  the  modified  SI  model 
[17].  In  deterministic  settings,  the  system  has  been  exploited 
to  model  single  and  coupled  patch  populations  [18],  as  well 
as  testing  vaccine  strategies  [14,19].  Assume  that  the  popu¬ 
lation  is  sufficiently  large  so  that  the  various  subgroups  are 
assumed  to  be  continuous.  The  population  dynamics  is  de¬ 
scribed  by  susceptible  S{t)\  exposed,  but  not  yet  infectious, 
E{f)-,  infective  I{t).  The  recovered  R{i)  class  in  the  model  can 
be  derived  from  model  results  since  S3-E+l+R=\  [17]. 

Seasonality  is  input  into  the  model  via  the  contact  rate, 
p{t),  so  we  let  p(t)=PQ{\  +  d  cos  Ittf),  where  0=S(J<1. 
Other  parameters  used  to  quantify  the  dynamics  are  a  scep- 
tible  input  rate  pi,  (which  includes  the  birth  rate,  as  well  as  a 
possible  fixed  vaccine  control),  the  mean  latent  period,  a~\ 
and  the  infectious  period,  y~K  The  full  deterministic  rate 
equations  are  given  by 

^=PL[i+h(t)]-msi~fis, 
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dm 

dt 


W  T  T 

—  =  aE-yI-/jJ, 


(3) 


/?W  =  l-[YW+£W+/«], 

where  hit)  is  a  small  perturbation  used  for  vaccination.  That 
is,  when  hit)  is  negative,  the  input  of  susceptibles  into  the 
system  is  reduced.  Since  it  will  be  designed  to  be  adaptive 
stochastic  control,  hit)  will  also  depend  on  the  state  vari¬ 
ables. 

For  realistic  childhood  disease  parameters  chosen  here, 
theoretical  [20]  and  munerical  analysis  [17]  show  that  for 
almost  all  cases,  the  infective  and  exposed  population  follow 
each  other  in  time  to  first  order,  leading  to  a  reduction  which 
describes  a  modified  SI  model  (MSI),  given  by 

hit)]  -  fiSit)  -  mmsit), 
at 

(4) 

^  )y3(0/(0S(t)  -  ima)lit). 

dt  \(j.+  y/ 

The  parameters  used  for  measles  data  [20]  are  given  by  fi 
=0.02,  a=  1  /0.0279,  y=  1  /O.Ol,  /3o=  1575,  and  5=0.095,  and 
are  fixed  throughout  the  paper.  Here,  the  parameter  hit)  is  a 
time-dependent  vaccine  control  whose  value  we  will  calcu¬ 
late  adaptively,  and  depends  on  the  phase  space  location  of 
(S(/),/(t)). 

Following  the  discretized  stochastic  model  in  Eq.  (1),  we 
strobe  the  system  with  period- 1  to  create  a  Poincare  map. 
Without  loss  of  generality,  we  define  a  discrete  stochastic 
model  for  the  purposes  of  this  paper  [21].  Using  a  discrete 
stochastic  map  approach  will  allow  us  to  make  careful  and 
accurate  interpretations  in  terms  of  the  (5'(/),/(r))  variables, 
as  well  as  to  examine  the  interaction  of  the  dynamics  and 
control  with  the  underlying  topology  of  the  system.  We  con¬ 
sider  the  rmcontrolled  stochastic  system  (/i=0)  as  a  two- 
dimensional  map  F  of  a  region  D  into  itself 

(5,7)(r+l)  =  Fl(5,/)(r)]+r7(t),  (5) 

where  97  is  a  two-dimensional  random  variable  having  a  nor¬ 
mal  distribution  given  by  vix)=e~^^^  '^^^^/(27t]|2||’^^),  with 
S=diag  (rF),  and  we  choose  the  standard  deviation  to  be  cr 
=0.035.  Since  the  two-dimensional  deterministic  system  has 
an  attractor  with  unequally-sized  components,  the  noise  am¬ 
plitude  is  scaled  so  that  it  is  defined  on  the  unit  square. 
Because  the  standard  deviation  is  based  on  the  rescaled  co¬ 
ordinates,  it  is  small  compared  to  the  attractor  size  and  is 
smaller  than  the  modulation  component  of  the  contact  rate  in 
Eq.  (4).  A  typical  time  series  of  the  7  component  is  shown  in 
Fig.  1.  Notice  the  frequent  aperiodic  bursts,  which  for  the 
chosen  parameters  of  the  deterministic  part  of  the  model,  Eq. 
(4),  would  not  occur  were  it  not  for  the  random  perturbations 
in  Eq.  (5);  the  deterministic  and  stochastic  parts  interact  in  a 
fundamental  way  to  create  complicated  oscillations  that  ei- 


time 


FIG.  I.  This  is  an  uncontrolled  time  series  of  the  fraction  of 
infectives  /  for  the  MSI  model  under  random  forcing  in  Eqs.  (4)  and 
(5).  The  parameters  are  given  in  the  text.  Inset:  Note  the  small  (SA) 
period  2  and  large  (LA)  period  3  amplitude  oseillations  of  the  un¬ 
derlying  bistable  deterministic  system.  No  chaos  is  present  when 
7jit)=Q,  and  the  system  only  exhibits  periodic  SA  or  LA 
oscillations. 

ther  phenomenon  could  not  create  on  their  own. 

Notice  that  in  the  absence  of  any  stochastic  fluctuations 
[57(1)  =  0]>  the  system  will  settle  down  to  one  of  two  periodic 
solutions.  The  two  stable  solutions  are  plotted  in  the  Fig.  1 
inset.  The  period  2  cycle  has  a  small  amplitude  (SA)  while 
the  period  3  cycle  is  of  large  amplitude  (LA).  However,  as 
seen  from  the  time  series  in  the  figure,  outbreaks,  which 
occur  due  to  stochastic  fluctuations,  may  have  enhanced  am¬ 
plitudes  by  almost  an  order  of  magnitude  over  the  period  3 
cycle. 

Although  the  system  is  stochastic,  its  dynamics  may  be 
quantified  in  terms  of  Lyapunov  exponents  by  spatial  inte¬ 
gration  against  the  invariant  density  [22].  For  the  parameters 
used  to  generate  the  time  series  in  Fig.  1,  we  compute  the 
Lyapunov  exponents,  and  find  them  to  be  kj  =0.1 63 8  and 
X2=-0.4853.  These  values,  together  with  the  evidence  of 
nearly  intersecting  stable  and  unstable  manifolds  [23],  indi¬ 
cate  a  completed  horseshoe  dynamics  imder  the  influence  of 
the  noise,  described  as  stochastic  chaos  [8,15].  However,  the 
completed  horseshoe  dynamics,  indicative  of  chaos  in  deter¬ 
ministic  systems,  is  a  geometric  way  of  thinking  about  the 
interaction  of  noise  and  the  underlying  manifold  stmcture  of 
the  deterministic  part.  The  chaotic-looking  dynamics  are  the 
result  of  mixing  two  stable  attractors,  while  sampling  un¬ 
stable  dynamics  between  them.  The  positive  Lyapunov  expo¬ 
nent  is  therefore  a  way  of  measuring  contributions  to  the 
stochastic  attractor  of  dynamics  tracking  near  imstable  mani¬ 
folds.  The  fraction  of  time  spent  near  the  unstable  manifolds, 
as  well  as  the  transition  probabilities  of  the  dynamics  switch¬ 
ing  from  small  to  large  amplitude  behavior  may  be  explained 
by  taking  a  dynamic  probabilistic  approach,  which  we  sketch 
briefly.  A  full  mathematical  description  is  given  in  [15]. 

ni.  DISCRETE  STOCHASTIC  DYNAMICS  AND 
TRANSITION  PROBABILITIES 

If  the  noise  is  continuous,  we  can  compute  the  evolution 
of  the  probability  density  using  a  Fokker-Planck  approach 
[24].  However,  since  the  approach  is  one  of  discrete  noise  as 
in  Eq.  (1),  we  evolve  the  densities  discretely  as  well.  That  is. 
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since  the  solution  to  the  periodically  driven  is  computed  ev¬ 
ery  period  to  form  the  discrete  map,  we  do  the  same  with  the 
density. 

We  assume  the  noise  comes  from  a  distribution,  v(x).  The 
evolution  of  an  initial  probability  density  function  (PDF), 
p:£)CR^— ►R,  is  defined  by  the  stochastic  Frobenius-Perron 
operator  [15]  Pf:L'(R^)-^f''(R^))  given  by 

^’f(pW)=  f  A.x-P'(y)]p(y)dy.  (6) 

J  D 

The  density  is  invariant  if  it  is  a  fixed  point  of  the  operator. 
This  approach  allows  an  approximation  of  the  probabilistic 
transitions  of  one  part  of  phase  space  to  another  [15]  as  well 
as  the  invariant  density  [25]. 

To  compute  the  transition  probabilities  from  one  region  of 
phase  space  to  another,  we  discretize  the  region  D  of  phase 
space.  Specifically,  we  assume  there  exists  a  cover  of  the 
region  D  by  disjoint  sets  5„ 

N 

D=Ufi,-.  (7) 

1=1 

Defining  the  set  of  characteristic  basis  functions. 


1,  X  e  Bi 
0,  X  ^  Bj 


(8) 


allows  one  to  generate  finite  dimensional  projections  of 
transport  by  computing  the  NX  N  matrix  entries  of  a  transi¬ 
tion  probability  matrix  [8,15]  given  by  the  equation 


^ij=  f  ^F(>/'i(x))</'j(x)dx. 
J  D 


Therefore,  Eq.  (9)  yields  the  probability  of  transporting  mass 
from  box  Bi  to  Bj. 

In  considering  the  problem  of  predicting  stochastic  out¬ 
breaks  in  the  MSI  model,  we  wish  to  compute  the  transition 
from  a  small  amplitude  (SA)  oscillation  to  a  large  amplitude 
(LA)  outbreak  in  a  time  series,  such  as  the  one  generated  in 
Fig.  1.  The  inset  shows  the  deterministic  periodic  orbits  of 
SA  and  LA,  although  noise  may  generate  much  larger  out¬ 
breaks  than  the  deterministic  LA  orbit.  Stochastic  perturba¬ 
tions  of  SA  in  the  inset  are  approximately  the  same  ampli¬ 
tude,  and  therefore  are  used  as  a  threshold  to  define  large 
outbreaks.  The  mass  flux  entries  generated  by  Eq.  (9)  can  be 
combined  with  the  invariant  density  to  generate  the  condi¬ 
tional  probability  of  transition  from  set  S,-  to  Bp  given  fi,-.  A 
representation  of  the  transition  probability  is  depicted  in  Fig. 
2.  Notice  that  the  most  active  transport  regions  lie  close  to  a 
stable  manifold  of  an  LA  orbit  (period  3  saddle)  in  the  un¬ 
derlying  deterministic  system.  This  stable  manifold  is  the 
deterministic  basin  boxmdary  which  separates  the  SA  (period 
2)  and  LA  (period  3)  regular  orbits  of  Eq.  (4),  and  the  col¬ 
oring  denotes  the  degree  and  location  where  this  pseudobar¬ 
rier  is  overcome  due  to  noise.  Notice  that  near  each  of  the 
basin  boundary  saddles  of  period  3,  transition  to  an  outbreak 
is  likely.  However,  the  highest  transition  region  is  not  near 
any  saddle.  Rather,  the  probability  of  an  outbreak  in  this 


region  is  solely  due  to  the  interaction  of  the  noise  and  the 
global  topology  of  the  underlying  deterministic  dynamics. 

IV.  ACTIVE  CONTROL  OF  STOCHASTIC 
OUTBREAKS 

For  deterministic  systems,  normal  methods  of  vaccine 
control  will  reduce  the  input  rate  of  susceptibles.  The  value 
of  h  is  usually  computed  so  that  at  equilibrium  (no  seasonal 
forcing,  or  5=0),  the  net  rate  of  production  of  infectives  in 
one  infectious  period  is  less  than  unity.  Under  these  condi¬ 
tions,  the  disease  will  die  out.  However,  control  of  small 
amplitude  oscillations  in  the  periodically  driven  case  can  be 
done,  but  the  disease  will  persist  [26]. 

In  the  stochastic  case  in  the  presence  of  periodic  drives, 
constant  controls  may  make  the  problem  worse.  In  Fig.  3,  we 
see  a  direct  comparison  of  constant  vaccine  control  and  no 
control.  Notice  Aat  although  the  mean  level  of  outbreaks 
appear  to  be  reduced,  the  large  fluctuations  are  greater  than 
without  control.  Therefore,  constant  vaccine  control,  al¬ 
though  sometimes  the  only  guide,  may  increase  the  size  of 
large  outbreaks.  Therefore,  it  is  natural  to  try  to  sample  and 
control  discretelywhen  considering  stochastic  outbreaks. 

Vaccine  activation  using  a  variable  h  depends  on  finding 
the  regions  where  an  outbreak  is  most  likely  upon  the  next 
iteration.  These  are  points  of  the  trajectory  generated  by  Eq. 
(5)  in  the  SA  basin  that  precede  iterates  in  the  LA  basin. 
Although  we  compute  conditional  outbreaks  from  the  spatial 
averages  using  the  transition  matrix,  this  is  verified  tempo¬ 
rally.  Using  an  uncontrolled  stochastic  time  series  of  50  000 
iterates,  and  checking  in  which  basin  (SA  or  LA)  each  iterate 
is  located,  we  show  in  Fig.  4  the  most  likely  preoutbreak 
regions.  In  comparison  to  Fig.  2,  the  spatial  average  predicts 
similar  transport  regions  of  high  conditional  probability  of 
the  SA-LA  transition. 

We  now  define  a  bull’s  eye  (BE)  region  to  be  an  open 
connected  neighborhood  having  high  probability  of  transi¬ 
tion  from  SA  to  LA  outbreaks.  The  BE  region,  for  a  chosen 
threshold,  is  clearly  shown  in  red  in  Fig.  2.  Distinguishing 
the  center  point  x^.,  the  BE  region  includes  a  neighborhood  of 
radius  e  that  has  a  probability  greater  than  a  given  threshold. 
Notice  that  this  is  not  the  only  region  in  which  transition 
occurs.  Monitoring  the  BE  region  alone,  therefore,  is  not 
sufficient  for  prediction  of  transitions  [27].  However,  it  can 
be  used  to  determine  other  regions  that  are  not  obvious  for 
transition  to  an  outbreak. 

We  can  use  the  BE  region  as  a  first  guess  to  monitor  the 
dynamics.  Let  xq  e  £  be  the  current  point  of  the  observed 
dynamics,  and  a  desired  target  point  in  the  transport  space 
close  to  the  image  of  E,  but  in  a  region  of  lower  transition 
probability.  The  relationship  between  the  current  point  in  the 
trajectory  xq  and  the  center  of  the  bull’s  eye  x^  is  XQ=Xc+y 
for  some  y.  To  move  the  image  of  Xo  closer  to  the  target  point 
xl,  we  activate  the  control  parameter  h  in  Eq.  (4).  By  Taylor 
expansion  about  F(X(,,/i)  when  h=0  and  ignoring  higher  or¬ 
der  terms,  we  solve 


[xi  -  Fjxc/x)  -  d^FiXc,fi)yYdJ^{x^,fi) 
\\d^F{x^,lx)f 


assuming  dpF{xi.,fi)  #  0.  This  control  strategy  is  designed  to 
target  a  desired  region  of  lower  probability,  given  the  iterates 
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FIG.  2.  (Color)  The  GTM  result  of  the  conditional  probability  of  transition  from  small  amplitudes  to  large  outbreaks  using  the  same 
parameters  as  in  Fig.  1.  The  highest  probability  regions  of  transport  (red)  point  to  a  bull’s  eye  monitoring  region  for  control.  Overlaid  are  the 
stable  and  unstable  manifolds  corresponding  to  the  underlying  deterministic  model. 


land  in  a  region  of  high  transition  probability. 

Now  we  apply  control  to  suppress  large  amplitude  out¬ 
breaks.  Focusing  on  points  in  a  neighborhood  of  the  BE  re¬ 
gion  has  the  disadvantage  that  the  values  of  I  are  already 
fairly  large.  Therefore,  we  use  the  detection  region  of  the 
neighborhood  around  the  (deterministic)  preimage  of  the 
bull’s  eye,  shown  as  an  ellipse  in  Fig.  4.  Using 

Eq.  (10),  the  image  of  the  ellipse,  //,,  is  found  to  be  the  figure 
eight  shown  in  Fig.  4.  We  targeted  a  region  in  4  which  is 
close  to  the  BE  region  but  has  a  very  low  transition  probabil¬ 
ity.  Our  techniques  successfully  steer  trajectories  away  from 
the  bull’s  eye  region  towards  SA  behavior  by  using  only 
vaccine  perturbations  that  control  the  flow  of  susceptibles 
about  some  mean  value. 

One  advantage  of  choosing  the  detection  region  to  be  the 
preimage  of  BE  is  for  relatively  low  values  for  the  number  of 
infected  individuals  (/),  a  prediction  can  be  made  about  the 
future  increase  and  steps  can  be  taken  to  avert  these  dynam¬ 
ics.  The  perturbations  represent  a  vaccination  program,  tak¬ 
ing  the  form  of  At„ew=At[l  +A(r)].  If  h  is  negative,  then  more 
vaccinations  are  required  to  reduce  the  rate  of  susceptible 
individuals  being  introduced  into  the  population.  An  example 


of  the  success  of  this  algorithm  is  shown  in  Fig.  5.  On  aver¬ 
age,  perturbations  are  applied  25-30  %  of  the  time.  Notice 
the  maximum  amplitude  in  comparison  to  the  uncontrolled 
dynamics  of  Fig.  1 .  For  this  example,  the  Lyapunov  expo¬ 
nents  are  X.| =0.0794  and  X.2~“0'3764,  where  the  maximum 
exponent  has  been  significantly  decreased. 

V,  DISCUSSION 

Stochastic  bursting  is  present  in  many  systems  that  are 
based  on  population  dynamic  modeling.  In  general,  when 
such  systems  are  subject  to  periodic  forcing,  there  exist  pa¬ 
rameter  regions  in  which  multiple  attractors  coexist.  Typi¬ 
cally,  one  of  these  attractors  arises  from  periodically  forced 
equilibrium,  and  therefore,  is  typically  of  small  amplitude. 
On  the  other  hand,  the  other  attractors  bifurcate  from  saddle 
node  orbits,  which  tend  to  be  of  larger  amplitude.  Such 
bistable  systems  can  have  a  simple  manifold  structure,  but 
when  considered  in  the  presence  of  stochastic  fluctuations, 
they  may  exhibit  complex  mixing  between  the  bistable  at¬ 
tractors,  coupled  with  complicated  looking  transients  be¬ 
tween  the  basins. 
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time 


FIG.  3.  (a)  An  uncontrolled  time  series  of  infective  fraction  as  a  function  of  time,  (b)  Constant  vaccine  control  to  reduce  the  rate  of  input 
of  susceptibles. 


By  using  the  PDF  flux,  we  are  able  to  distinguish  regions 
in  the  small  amplitude  basin  that  are  quite  sensitive  to  sto¬ 
chastic  effects.  We  use  this  information  in  a  control  algo¬ 
rithm  to  prevent  bursting  dynamics  (that  is,  to  control  sto¬ 
chastic  chaos).  It  monitors  this  sensitive  region  and  adjusts 
one  physically  relevant  parameter  to  keep  trajectories  in  the 
SA  basin.  This  idea  of  monitoring  a  loss  region  has  been 


-3.0  -2.7  -2.4 

FIG.  4.  Temporal  average  of  those  iterates  leading  to  outbreaks 
in  the  next  iterate  using  the  same  parameters  as  in  Fig.  1.  Notice  the 
agreement  with  the  spatial  average  in  Fig.  3.  The  ellipse  bounds  the 
detection  region.  The  figure  eight  curve  is  the  image  of  the  ellipse 
with  controlled  targeting. 


used  in  other  chaos  control  schemes  that  are  deterministic 
(i.e.,  [28,29]).  To  our  knowledge,  we  are  not  aware  of  any 
stochastic  chaos  control  methods  that  account  specifically  for 
the  emergent  effects  of  stochastic  perturbations. 


FIG.  5.  Stochastic  control  to  suppress  large  outbreaks  in  the 
MSI  model,  (a)  Infectives  with  suppressed  outbreaks  due  to  control 
in  the  influx  of  infectives.  (b)  Perturbations  h  to  the  susceptible 
input  rate  /j,  in  Eq.  (4). 
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FIG.  6.  A  plot  showing  optimized  predictability.  Plotted  are  the 
false  alarm  rates  versus  missed  detections. 

One  concern  with  a  probabilistic  detection  scheme  is  that 
it  is  dependent  on  the  choice  of  the  monitoring  region  used 
for  transition  to  an  outbreak.  Two  issues  with  taking  an  ac¬ 
tual  time  series  and  using  the  monitoring  scheme  above  is 
that  it  may  miss  an  outbreak  that  is  there  (missed  detection), 
or  it  may  predict  an  outbreak  that  does  not  occur.  These 
statistics  depend  heavily  on  the  size  of  the  monitoring  one 
uses.  To  see  this  in  Fig.  6,  we  change  the  radius  around  of  the 
center  of  the  bull’s  eye  and  the  radius  around  its  preimage. 
Each  dot  plotted  in  Fig.  6  is  for  a  different  radius.  The  small¬ 
est  radii  are  represented  by  the  data  points  on  the  right.  As 
we  increase  the  radii,  the  data  points  move  along  the  curve  to 
the  left.  The  false  alarms  are  those  outbreaks  predicted  by  the 
bull’s  eye,  but  do  not  occur.  The  missed  detection  are  the 
bursts  that  occur  but  are  not  predicted  by  the  maximum  flux 
hypothesis.  It  is  the  percentage  not  detected. 

The  choice  we  made  for  the  detection  region  has  solely 
been  guided  by  time  series  observations  and  PDF  flux  pre¬ 
dictions.  It  has  not  been  optimized  for  the  minimum  number 
or  size  of  perturbations.  Because  of  the  stochastic  perturba¬ 
tions  added  to  the  system,  the  control  measures  will  not  “trail 
off”  as  in  targeting  unstable  periodic  orbits. 

One  very  interesting  aspect  of  the  control  perturbation  is 
that  the  value  of  h  is  sometimes  positive.  This  counterintui¬ 
tive  result  may  be  explained  since  an  increase  in  S  is  used  to 


trigger  an  earlier,  but  smaller,  outbreak.  To  xmderstand  this, 
we  consider  the  MSI  model,  but  transformed  and  scaled,  so 
that  the  steady  state  equilibrium  in  the  absence  of  forcing  is 
now  at  the  origin,  and  we  examine  the  conservative  system 
in  the  absence  of  damping  as  well  [30], 


x'{t)  =  -vy, 

(11) 

/(/)  =  wc(l  +y). 

In  Eq.  (II),  X  is  a  scaled  susceptible,  y  is  a  scaled  infec¬ 
tive,  the  equilibrium  is  at  the  origin,  and  the  frequency  v  is  a 
function  of  the  epidemiological  parameters.  Notice  that  since 
the  population  is  assumed  to  be  constant,  inthe  absence  of 
any  infectives  [><=-1  in  Eq.  (11)],  the  fraction  of  suscep- 
tibles  slowly  increases.  In  addition,  all  oscillatory  solutions 
must  lie  on  level  curves  to  the  Lyapunov  function:  V(x,y) 
=x2+2y-2  1n(y+l). 

Now  suppose  we  have  a  small  amount  of  infectives  im¬ 
posed  by  a  strong  level  of  vaccine.  Then  the  infectives  will 
stay  small  for  a  long  period  of  time,  until  enough  suscep- 
tibles  grow  to  cause  an  outbreak  of  very  large  amplitudes  by 
coming  in  contact  with  a  few  infectives  [31].  That  is,  an 
outbreak  will  not  occur  unless  the  susceptibles  reach  a  criti¬ 
cal  level  in  a  long  time  scale  while  in  the  presence  of  a  small 
fraction  of  infectives.  To  be  specific,  suppose  y=— l+ce, 
where  c>0  is  constant.  Theny'=xce  If  x<0,  then  the  in¬ 
fectives  decrease  further,  implying  a  much  larger  outbreak  at 
a  later  time.  Therefore,  if  one  increases  the  infectives,  the 
system  fires  sooner,  with  a  smaller  outbreak,  since  the  infec¬ 
tives  are  pushed  further  away  from  the  invariant  line  y=-l. 
When  the  control  h  is  adjusted  so  that  it  is  positive,  the  effect 
is  to  cause  an  increase  in  the  rate  of  infectives,  thus  reducing 
the  size  of  the  outbreak. 

Finally,  although  the  vaccine  control  fluctuations  do  not 
decrease  the  mean  incidence  levels  of  infection,  the  control 
may  be  combined  with  tracking  methods  for  epidemic  con¬ 
trol  [26]  to  reduce  the  mean  reproductive  rate  of  infection 
below  threshold  to  kill  off  the  disease  without  causing  un¬ 
wanted  outbreaks  during  vaccination. 
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We  analyze  the  effects  of  stochastic  perturbations  in  a  physical  example  occurring  as  a  higher-dimensional 
dynamical  system.  The  physical  model  is  that  of  a  class-5  laser,  which  is  perturbed  stochastically  with  finite 
noise.  The  effect  of  the  noise  perturbations  on  the  dynamics  is  shown  to  change  the  qualitative  nature  of  the 
dynamics  experimentally  from  a  stochastic  periodic  attractor  to  one  of  chaoslike  behavior,  or  noise-induced 
chaos.  To  analyze  the  qualitative  change,  we  apply  the  technique  of  the  stochastic  Frobenius-Perron  operator 
[L.  Billings  et  al,  Phys.  Rev.  Lett.  88,  234101  (2002)]  to  a  model  of  the  experimental  system.  Our  main  result 
is  the  identification  of  a  global  mechanism  to  induce  chaoslike  behavior  by  adding  stochastic  perturbations  in 
a  realistic  model  system  of  an  optics  experiment.  In  quantifying  the  stochastic  bifurcation,  we  have  computed 
a  transition  matrix  describing  the  probability  of  transport  from  one  region  of  phase  space  to  another,  which 
approximates  the  stochastic  Frobenius-Perron  operator.  This  mechanism  depends  on  both  the  standard  devia¬ 
tion  of  the  noise  and  the  global  topology  of  the  system.  Our  result  pinpoints  regions  of  stochastic  transport 
whereby  topological  deterministic  dynamics  subjected  to  sufficient  noise  results  in  noise-induced  chaos  in  both 
theory  and  experiment. 
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I.  INTRODUCTION 

Noise-induced  escape,  which  appears  as  a  form  of  bifur¬ 
cation  in  dynamical  systems,  is  now  documented  in  many 
areas  of  science  and  engineering  [1].  It  arises  in  stochastic 
processes,  which  we  consider  to  be  a  composition  of  deter¬ 
ministic  and  time-dependent  noisy  systems.  Detecting  chaos 
in  noisy  systems  is  still  an  issue  of  debate.  Efforts  have  been 
made  to  carry  over  operational  definitions  of  deterministic 
chaos  to  stochastic  systems,  such  as  proving  the  existence  of 
a  positive  Lyapunov  exponent  [2]  and  exploring  the  interac¬ 
tion  of  noise  and  a  global  bifurcation  based  on  underlying 
unstable  structures,  such  as  a  chaotic  saddle  [3].  Many  of  the 
underlying  deterministic  systems  in  these  examples  have  pa¬ 
rameter  regimes  in  which  multiple  attractors  give  rise  to 
noise-induced  escape  from  one  attractor  to  another.  Such  sys¬ 
tems  may  be  analyzed  globally  using  the  Hamiltonian  theory 
of  large  fluctuations,  or  considering  escape  from  attracting 
potential  wells  along  most  probable  exit  paths  [4]  using  the 
theory  of  quasipotentials  [5,6]  or  a  variational  formulation  of 
optimal  escape  paths  [7].  It  is  well  known  that  noise  can 
excite  unstable  chaotic  structures  while  destroying  regular 
periodic  dynamics,  but  most  studies  consider  noise-induced 
chaos  occurring  near  a  bifurcation,  such  as  a  saddle-node 


point  or  a  crisis  of  chaotic  attractors  that  leave  a  chaotic 
saddle  present. 

In  driven  deterministic  systems,  the  existence  of  chaotic 
invariant  sets,  such  as  chaotic  saddles,  can  be  proven  by 
examining  the  topology  of  intersecting  manifolds  [8].  As  an 
example,  we  cite  the  Melnikov  method  [9].  Although  it  has 
been  extended  to  stochastic  systems  [10],  it  is  limited  in 
application  since  it  is  a  bifurcation  result  that  is  perturbed 
from  a  global  homoclinic  or  heteroclinic  connection  in  a  con¬ 
servative  system.  Therefore,  in  many  cases,  one  must  rely  on 
algorithmic  methods  for  the  numerical  computation  of  un¬ 
stable  objects  and  their  manifolds  [11-13],  with  the  hope  that 
one  may  extract  transverse  intersections.  We  also  note  that  in 
contrast  to  the  hypothesis  that  noise-induced  chaos  is  caused 
by  a  chaotic  saddle  excitation,  a  recent  result  shows  that  only 
partially  formed  manifold  intersections  (in  which  no  chaotic 
saddle  exists)  may  also  be  found  to  have  a  positive  Lyapimov 
exponent  [14]. 

In  this  paper,  we  compare  a  bifurcation  observed  in  a 
nonequilibrium  stochastic  class-5  laser  experiment  to  a  cor¬ 
responding  model  of  the  system.  We  include  experimental 
results,  as  well  as  the  theoretical  explanation  of  the  observa¬ 
tions.  In  particular,  experiments  support  the  claim  that  add¬ 
ing  larger  stochastic  perturbations  to  the  system  results  in 
qualitatively  different  dynamics.  Using  the  model,  we  pro- 
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FIG.  1.  Experimental  apparatus  used 
to  perform  our  measurements.  AOM, 
acousto-optic  modulator;  BS,  beam  split¬ 
ter;  PM,  power  meter;  PD,  photodetec¬ 
tor;  TG,  trigger  generator;  ADC,  analog- 
to-digital  converter;  and  PC,  personal 
computer. 


vide  evidence  that  this  is  an  example  of  a  bifurcation  to 
noise-induced  chaos  by  explicitly  computing  the  probability 
transport  due  to  noise.  In  this  way,  the  interaction  of  noise 
and  the  underlying  topology  is  identified  in  the  emergent 
dynamics.  We  present  analytic  methods  that  specifically 
carry  out  the  task  of  constructing  the  invariant  density  and 
transition  probabilities  in  a  rigorous  manner  to  address  the 
problem  of  this  P-type  stochastic  bifurcation,  as  defined  in 
[15].  Since  it  is  a  global  approach,  it  is  an  alternative  to  using 
the  Hamiltonian  theory  of  large  fluctuations,  as  described  in 
[7,16]  for  autonomous  systems,  and  in  [17]  for  periodic  sys¬ 
tems.  New  tools  were  developed  that  are  based  on  discrete 
approximations  to  the  Frobenius-Perron  operator  with  addi¬ 
tive  noise,  defined  as  the  stochastic  Frobenius-Perron  opera¬ 
tor  (SFPO)  [18,19]. 

Using  the  SFPO,  we  identify  the  active  regions  of  sto¬ 
chastic  transport,  or  probability  transitions,  in  the  model.  The 
advantage  of  this  method  is  that  we  can  find  the  probability 
density  function  (PDF)  and  maximal  transport  across  bound¬ 
aries  in  the  absence  of  a  priori  knowledge  of  manifold  struc¬ 
tures  and  without  time  averaging.  From  the  SFPO  method, 
since  one  can  directly  compute  the  invariant  density,  spa¬ 
tially  averaged  Lyapunov  spectra  may  be  computed  if  the 
linear  variation  along  an  orbit  is  known.  For  stochastic  sys¬ 
tems  that  are  sufficiently  ergodic,  spatial  and  temporal  aver¬ 
ages  of  the  Lyapunov  spectra  are  equal,  and  therefore,  a  posi¬ 
tive  Lyapunov  exponent  averaged  spatially  is  a  possible 
indicator  of  stochastic,  or  noise-induced,  chaos  [15]. 

To  contrast  our  work  from  previous  theories,  we  note  that 
the  bifurcation  is  far  from  parameters  that  would  lead  to  a 
natural  bifurcation  to  chaos,  and  large  noise  levels  are  in¬ 
cluded.  Many  studies  in  this  field  have  relied  on  examining 
small  noise  limits,  such  as  quasipotential  theory  [6]  and  op¬ 
timal  path  theory  [4],  although  this  work  has  more  recently 
been  extended  to  the  regime  of  finite  noise  intensity  [20]. 
Underlying  unstable  fractal  structures  and  noise-induced  ba¬ 
sin  escape  times  have  also  been  examined  from  quasipoten¬ 
tial  theory  [6]  for  simple  maps.  The  basin  boundary  in  the 
system  we  study  is  a  simple  stmcture;  i.e.,  it  is  nonffactal 
due  to  the  lack  of  intersecting  stable  and  unstable  manifolds. 
In  fact,  only  the  forward  crossings  of  a  heteroclinic  tangle 
could  be  identified,  and  no  nonattracting  chaotic  sets  are 
found  to  exist.  The  maximum  Lyapunov  exponent  was  cal¬ 
culated  to  increase  smoothly  through  zero  at  the  transition. 
Although  both  smooth  and  discontinuous  onset  are  attributed 
to  noise-induced  chaos,  the  transition,  which  resembles  a 
noise-induced  attractor  explosion  described  in  [6],  is  not  ob¬ 


served.  In  fact,  the  transition  is  smooth  statistically,  evi¬ 
denced  by  the  smooth  transition  of  a  Lyapunov  exponent 
through  zero,  which  may  be  due  to  the  noise-induced  un¬ 
stable  dimension  variability  [21].  There  is  also  a  resemblance 
to  noise-induced  switching  between  multiple  attractors  as  de¬ 
scribed  in  [3],  but  evidence  provided  by  the  probability  den¬ 
sity  fimction  supports  the  fact  that  trajectories  spend  as  much 
time  (if  not  more)  near  the  partially  formed  heteroclinic 
tangle  as  the  two  periodic  attractors.  We  also  note  that  as 
previously  reported  in  [6]  and  [3],  both  explosions  and  at¬ 
tractor  switching  are  facilitated  by  fractal  basin  boimdaries 
and  nonattracting  chaotic  sets. 

The  layout  of  the  paper  is  as  follows.  In  Sec.  II,  we  de¬ 
scribe  the  experimental  setup  of  a  nonequilibrium  stochastic 
class-5  laser.  We  illustrate  the  effects  of  noise  on  the  dynam¬ 
ics  of  the  intensity  and  show  how  the  stmcture  of  the  attrac¬ 
tor  changes.  In  Sec.  Ill,  we  briefly  review  the  laser  model  of 
the  experiment  in  a  reduced  form  and  show  that  it  captures 
many  of  the  features  of  the  experiment.  Section  IV  illustrates 
the  effect  of  noise  on  the  laser  model  and  specifically  shows 
how  the  maximal  (or  top)  Lyapunov  exponent  depends 
smoothly  on  the  standard  deviation  as  it  transitions  from  sto¬ 
chastic  periodic  behavior  to  stochastic  chaos.  The  global 
stmcture  of  the  xmderlying  topology  and  transport  results  are 
presented  in  Sec.  V,  and  the  discussion  is  presented  in  Sec. 
VI. 

II.  AN  ACOUSTICALLY  OPTICAL  MODULATED  LASER 
EXPERIMENT  WITH  NOISE 

To  examine  the  effects  of  external  noise  in  an  experiment, 
we  consider  an  acoustically  optical  modulated  laser  system. 
The  experimental  apparatus  is  shown  in  Fig.  1 .  It  consists  of 
a  single-mode  CO2  laser  with  an  intracavity  acousto-optic 
modulator  allowing  modulation  of  the  cavity  losses.  The  op¬ 
tical  cavity  is  1.30  m  long  and  the  total  transmission  coeffi¬ 
cient  r  is  0.10  for  a  single  pass.  The  intensity  decay  rate  k{i) 
can  be  expressed  as  follows: 

k{t)  =k{\+ a  sin2{5o[  1  +/(t)]}) ,  ( 1 ) 

where  k=cTIL,  c  is  the  speed  of  light  in  a  vacuiun,  L  is  the 
cavity  length,  q'=(1-27)/27’,  Bq  is  a  bias,  andy[t)  is  the 
modulation  signal, 

f{t)  =  13  sin(2irvf)  +  77(1) ,  (2) 

with  v=100  kHz  and  the  modulation  amplitude  The  ran¬ 
dom  variable  rj  is  considered  to  be  normally  distributed  with 
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FIG.  2.  Contour  plots  of  the  embedded  intensity  data  (arbitrary 
units)  under  perturbations  of  t).  Darker  shades  indicate  regions  vis¬ 
ited  with  higher  frequency.  Small  perturbations  are  used  in  the  left 
graph,  which  result  in  a  noisy  periodic  orbit.  Larger  perturbations 
are  used  in  the  right  graph.  Notice  how  the  emergent  dynamics  are 
fundamentally  different  from  the  smaller  case. 

mean  zero  and  standard  deviation  cr.  The  noisy  signal  f{{)  is 
provided  by  an  arbitrary  waveform  generator  (Tektronix 
Mod.  AWG420),  which  generates  both  the  sinusoidal  signal 
and  the  random  variable  rj  using  an  independent  internal 
Gaussian  noise  generator.  Specifically,  the  noise  is  added  pe¬ 
riodically  with  the  period  of  drive;  i.e., 

=  /i  =  l,2, ...,  (3) 

and  S  is  the  Dirac  delta  function,  and  7]„  is  now  a  discrete 
random  variable. 

It  is  known  that  by  increasing  the  amplitude  modulation, 
the  system  undergoes  a  sequence  of  subharmonic  bifurca¬ 
tions  leading  to  chaos  when  17= 0  [22].  However,  when  noise 
is  added  to  the  system  through  the  driver,  the  resulting  dy¬ 
namics  is  highly  dependent  on  the  noise  amplitude.  In  Fig.  2, 
we  see  two  examples  of  the  output  of  the  intensity  plotted  as 
a  contour  map  of  the  embedded  data  for  two  values  of  the 
noise  strength  at  the  same  value  of  the  modulation  amplitude 
y8=0.360.  Note  that  darker  shades  indicate  regions  visited 
with  higher  frequency,  and  x„  is  the  local  maximum  of  the 
measured  intensity.  The  left  panel  shows  the  case  where 
small  noise  results  in  a  two-piece  attractor.  The  deterministic 
attractor  at  this  parameter  value  is  periodic,  located  at  the 
dark  regions  in  the  middle  of  the  pieces.  We  describe  this 
behavior  as  a  noisy  periodic  attractor  since  its  power  spec¬ 


trum  is  slightly  broadened  around  the  corresponding  fre¬ 
quency.  In  the  right  panel,  considerably  more  noise  causes  a 
qualitative  change  in  the  attractor.  The  periodic  orbit  is  still 
located  in  the  darkest  regions  of  the  graphs,  but  notice  how 
there  is  significant  sampling  to  other  parts  of  the  phase  space 
not  previously  sampled  at  lower  noise  amplitudes.  To  rmder- 
stand  the  mechanism  resulting  in  the  stochastic  bifurcation, 
we  consider  an  accurate  model  of  the  experimental  process. 

HI.  AN  AOM  LASER  MODEL 

In  [23],  a  multifrequency  phase  control  on  a  two-level, 
two-dimensional  CO2  laser  model  produced  both  experimen¬ 
tal  and  numerical  evidence  that  it  was  able  to  preserve  peri¬ 
odic  behavior  within  a  chaotic  window  as  well  as  to  reexcite 
chaotic  behavior  when  it  is  destroyed  by  a  crisis.  In  the 
model  used,  only  intensity  and  population  inversion  were 
considered.  To  retain  fidelity  between  theory  and  experiment, 
a  more  realistic  four-level  model  of  a  CO2  laser,  which  in¬ 
corporated  intensity,  two  resonant  population  levels,  and  two 
coupled  rotational  manifolds,  was  used  in  [24].  Analysis 
showed  that  an  approximate  reduction  to  three  state  variables 
could  be  made  by  examining  differences  in  the  resonant  and 
rotational  population  levels  while  still  retaining  many  of  the 
global  features  of  the  bifurcations.  Therefore,  we  begin  our 
study  of  the  scaled  three-dimensional  model  in  a  stochastic 
version,  where  noise  is  added  to  the  intensity  equation.  The 
variables  have  already  been  scaled  to  be  dimensionless  [24]. 
The  driven  three-dimensional  system  has  the  advantages  that 
(i)  it  is  higher  dimensional  than  other  models,  and  (ii)  when 
sampled  discretely  at  the  drive  frequency,  its  phase  space  can 
be  visualized  in  three  dimensions.  The  model  equations  are 
given  by 

y\  =  *o(y2  -  1  -  a  sin^{5[l  +/(/)]}), 
y2  =  -  71T2  -  2A:o^>>'2  +73  + 

73  =  “7273+272+2^.  (4) 

and 

/(f)=.4  sin(wf)+ j/t),  (5) 

where  77(f)  is  discretely  modeled  as  in  Eq.  (3)  with  period 
2'itI  (o,  yi  is  the  natural  logarithm  of  the  intensity,  ^2  is  the 
main  population  difference,  and  is  the  difference  in  rota¬ 
tional  levels.  The  fixed  parameters  are  kQ=2)2.91,  a;=4,  B 
=0.21,  0=0.897  597  9,  ri  =  10.0643,  P=0.082,  72=1  0643, 
z=10,  and  we  vary  A. 

We  now  describe  the  topology  of  Eqs.  (4)  and  (5)  without 
stochastic  perturbations,  i.e.,  7^t)  =  Q.  As  shown  in  the  bifur¬ 
cation  diagram  in  Fig.  3,  periodic  orbits  are  represented  as  a 
function  of  A.  As  A  is  increased,  a  period-one  attractor  pro¬ 
ceeds  through  a  period-doubling  bifiircation.  Several  saddle- 
node  bifurcations  for  varying  periodic  orbits  also  occur, 
which  will  play  a  role  when  noise  is  turned  on.  We  show  the 
first  saddle  node,  which  is  of  period  three  in  the  figure. 
Therefore  at  vf  =0.214,  there  exists  an  interval  of  bistability 
formed  by  period-four  and  period-three  attractors.  Asso- 
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FIG.  3.  (Color  online)  The  bifurcation  diagram  for  the  laser 
model  as  a  function  of  the  forcing  amplitude.  Plotted  are  branches 
of  both  stable  and  unstable  periodic  orbits.  The  y  axis  is  scaled 
intensity.  The  parameters  are  given  in  the  text. 

ciated  with  the  period-four  attractor  is  an  unstable  period- 
two  saddle  orbit  (which  is  a  flip  saddle)  and  an  unstable 
period-one  flip  saddle  orbit  from  the  period-doubling  bifur¬ 
cation.  The  period-three  attractor  has  an  associated  unstable 
period-three  regular  saddle  orbit  arising  from  a  saddle-node 
bifurcation.  We  hypothesize  that  the  multi-instability  in  this 
system  when  Ae.1^  has  the  topological  structure  needed  to 
induce  chaoslike  behavior  with  additive  stochastic  perturba¬ 
tions.  Since  the  bifurcation  diagram  contains  only  stability 
and  amplitude  information,  we  explore  the  phase  space 
through  numerical  simulation.  (Rigorous  analysis  of  the  on¬ 
set  of  the  saddle-node  bifurcation,  which  leads  to  bistable 
regions,  is  similar  to  that  done  in  [25]  and  will  be  presented 
elsewhere.) 

IV.  STOCHASTIC  DYNAMIC  SIMULATIONS 

In  keeping  with  the  experimental  setup,  we  model  the 
stochastic  system  as  a  discrete  dynamical  system.  Since  the 
experimental  system  was  forced  periodically  with  discrete 
noise  using  Eq.  (3),  we  can  add  the  perturbations  at  the  same 
period  as  that  of  the  drive  given  by  Eq.  (5).  Consider  the 
periodic  sampling  as  discrete  time  events  of  a  deterministic 
system.  We  add  the  perturbations  to  initial  conditions,  similar 
to  adding  noise  to  a  discrete  map.  In  general,  we  consider 
stochastically  perturbing  a  function  F  with  additive  noise; 
F:R^— +R^,  xi->F(x)+ where  tj  is  an  identically  indepen¬ 
dently  distributed  random  variable  with  normal  distribution 
and  mean  =0  applied  once  each  iteration.  Since  we  are  most 
interested  in  the  situation  where  small  noise  amplitude  can 
have  major  global  consequences,  we  focus  on  the  case  where 
the  random  part  r)  is  assumed  to  be  independent  of  state  x 
and  relatively  small,  so  that  the  deterministic  part  F  has  pri¬ 
mary  influence.  We  add  the  perturbation  to  each  component 
independently  and  set  the  standard  deviation  2 
=diag(cr,  ,a-2,cr3)  as  a  parameter.  This  standard  deviation  is 
relative  to  the  normalized  scaling  of  the  almost-zorapact 
space  we  consider.  That  is,  each  (t„  is  scaled  as  if  the  phase 


space  is  a  unit  box  in  three  dimensions.  With  no  noise,  the 
only  observable  behavior  is  asymptotically  periodic  trajecto¬ 
ries  converging  to  the  period-three  or  period-four  orbits.  By 
adding  noise  with  increasing  standard  deviation,  a  random 
trajectory  changes  from  a  noisy  periodic  orbit  to  chaoslike 
deterministic  behavior,  visiting  the  two  periodic  orbit  basins. 

We  remark  that  although  in  the  original  imscaled  model, 
noise  is  added  multiplicatively,  it  is  approximately  equivalent 
to  adding  noise  additively  in  the  scaled  model  from  Eq.  (4). 
This  is  due  to  the  fact  that  the  intensity  of  the  original  model 
is  represented  by  the  logarithm  of  the  intensity  in  Eq.  (4) 
[26].  The  noise  source  appears  as  a  term  of  the  form 
sin^{R[l+y4  sin(ri)0  +  ’7(0]}-  Taking  a  Taylor  series  expan¬ 
sion  with  respect  to  r)  yields  a  noise  term  on  the  order  of  rj, 
which  is  independent  of  the  state  variables.  Since  the  model 
is  based  on  the  natural  logarithm  of  the  intensity,  a  good 
approximation  to  the  noise  source  is  that  the  intensity  equa¬ 
tion  has  an  additive  noise  term. 

In  quantifying  underlying  complex  determinism  in  sto¬ 
chastic  systems,  it  is  inherently  difficult  to  draw  a  clear  line 
to  distinguish  between  complex  oscillations  due  to  signifi¬ 
cant  contributions  from  deterministic  parts  influenced  by 
noise  and  a  large  noise  amplitude  effect  wherein  complex 
oscillations  are  primarily  due  to  random  Brownian  diffusion. 
One  necessary,  but  not  sufficient,  condition  for  the  existence 
of  chaos  is  the  calculation  of  positive  Lyapunov  exponents. 
Lyapunov  exponents  measure  the  average  rate  of  separation 
of  neighboring  initial  points.  Because  we  are  adding  pertur¬ 
bations  to  this  system  discretely,  we  can  find  a  finite-time 
numerical  approximation  for  the  Lyapunov  exponents  of  the 
map  using  the  linear  variational  equations  of  the  original 
system  on  the  Poincare  section.  A  positive  Lyapunov  expo¬ 
nent  can  identify  chaotic  behavior,  but  diffusion  can  yield  a 
positive  Lyapunov  exponent  as  well  [27].  Since  chaos  is  also 
associated  with  the  underlying  topology  of  the  manifolds  of 
the  dynamical  system,  we  examine  the  unstable  stmctures  in 
the  deterministic  model  and  observe  how  they  interact  with 
the  stochastic  source  terms.  Specifically,  we  would  like  to 
identify  the  stmctures  in  the  original  phase  space  that  noise 
can  excite.  For  example,  if  noise  causes  a  trajectory  to  visit  a 
chaotic  saddle,  then  there  should  be  locally  unstable  contri¬ 
butions  to  the  Lyapunov  spectrum.  If  enough  of  the  xmstable 
contributions  are  sampled,  then  the  topology  imderlying  the 
chaotic  saddle  will  be  reflected  in  an  increasing  maximum 
exponent. 

In  modeling  the  experiment,  we  consider  additive  stochas¬ 
tic  perturbations  to  the  first  component,  setting  <r2=cr3=0. 
The  phase  space  projection  of  the  attractor  changes  qualita¬ 
tively  as  a  function  of  the  standard  deviation,  as  we  saw 
earlier.  However,  in  Fig.  4,  we  show  how  the  experiment  and 
model  both  appear  to  change  smoothly  as  the  standard  de¬ 
viation  increased.  This  is  reflected  in  the  time-averaged 
Lyapunov  exponent  computations.  That  is,  as  we  increase  cr\ 
away  from  zero,  the  Lyapunov  exponent  increases  and  has  a 
smooth  transition  from  negative  to  positive  values,  as  shown 
in  Fig.  5.  The  crossing  is  near  o-j =0.064.  As  an  example,  we 
graph  two  trajectories  to  show  the  emergent  dynamics  in  the 
three-dimensional  phase  space  in  Fig.  6.  Setting  (ri=0.04, 
the  largest  Lyapunov  exponent  is  negative,  predicting  noisy 
periodic  behavior,  as  seen  near  the  period-four  orbit.  Setting 
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FIG.  4.  The  left  graphs  show  intensity  data 
(arbitrary  units)  for  increasing  perturbations  of  rj. 
The  right  graphs  show  similar  results  for  the 
model  given  by  Eq.  (4).  Plotted  are  successive 
local  maxima  of  the  intensity  values.  Notice  that 
in  both  cases,  the  attractors  go  from  a  stochasti¬ 
cally  perturbed  period-four  cycle,  through  a  basin 
hopping  attractor,  and  then  to  bursting  among 
several  basins  of  attraction  from  the  deterministic 
case. 


(Ji=0.16,  the  largest  Lyapunov  exponent  is  positive,  predict¬ 
ing  chaoslike  behavior. 

More  detail  about  the  dynamics  can  be  obtained  by  cal¬ 
culating  the  bursting  statistics  as  a  function  of  the  standard 
deviation  of  the  noise.  We  approximate  the  burst  rate  by 
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FIG.  5.  The  largest  Lyapunov  exponent  as  a  function  of  the 
standard  deviation  of  the  noise.  The  transition  from  negative  to 
positive  values  is  smooth,  as  predicted  in  [21]. 


setting  a  threshold  for  at  0.009.  This  value  was  deter¬ 
mined  by  monitoring  a  trajectory  with  no  noise.  For  each 
standard  deviation  value,  we  count  the  number  of  points  in  a 
random  trajectory  above  the  threshold  and  divide  by  the  total 
number  of  points.  See  Fig.  7  for  the  results  using  trajectories 
700  000  points  long.  Notice  how  bursting  occurs  for 
cr>  0.12.  This  value  is  different  from  the  bifurcation  value 
predicted  by  the  Lyapunov  exponent.  Therefore,  we  will  in¬ 
vestigate  the  stochastic  dynamical  system  as  the  noise  pa¬ 
rameter  (T  is  varied. 

Experimentally,  we  observe  and  show  in  Fig.  6  the 
changes  that  occur  as  cr  varies.  There  exists  a  two-piece 
noisy  period- four  attractor  for  o' <0.064.  Then,  the  two 
pieces  join  into  one  attractor  for  0.064<o-<0.12,  which  is 
reflected  by  a  positive  Lyapunov  exponent.  Then,  the  trajec¬ 
tories  start  to  burst  and  visit  the  period  three.  This  statistic  is 
not  noticeable  until  a=0.\2.  The  amount  of  bursting  is  re¬ 
flected  in  the  burst  rate.  For  0.064<o'<0.12,  the  noise  pro¬ 
vides  the  transport  for  the  trajectories  to  visit  the  stable 
period-four  orbit,  the  unstable  period-two  orbit,  and  the  un¬ 
stable  period-one  orbit.  But  for  o'>0.12,  the  trajectory  visits 
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FIG.  6.  The  left  graph  shows  noisy  periodic  behavior  generated 
by  the  system  when  the  standard  deviation  of  the  noise  is  cri 
=0.04.  This  is  the  behavior  predicted  by  a  negative  Lyapunov  ex¬ 
ponent.  The  right  graph  shows  chaoslike  behavior  generated  by  the 
system  when  (ri=0.16.  This  is  the  behavior  predicted  by  a  positive 
Lyapunov  exponent  in  Fig.  5. 


both  the  period-three  and  period-four  attractors,  and  the 
manifolds  in  between  them.  See  the  right  graph  in  Fig.  6  as 
an  example.  We  will  now  explore  these  dynamical  changes 
by  a  transfer  operator-based  analysis,  and  we  will  compare 
the  results  to  the  topology  of  the  stable  and  unstable  mani¬ 
folds  of  the  corresponding  deterministic  system  and  interpret 
the  influence  of  the  added  noise. 
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FIG.  7.  The  burst  rate  as  a  function  of  the  standard  deviation  of 
the  noise.  Notice  that  the  bursting  increases  from  zero  near  a 
=0.12. 


V.  PHASE-SPACE  ANALYSIS  OF  STOCHASTIC 
DYNAMICS 

Understanding  the  interaction  between  noise  and  the  de¬ 
terministic  topology  requires  that  we  examine  the  structure 
of  the  stable  and  unstable  manifolds  of  the  relevant  saddles 
in  the  prechaotic  regime.  Locating  stable  and  unstable  mani¬ 
folds  can  be  done  in  several  ways  [11,12].  We  use  the  box 
algorithm  from  [13]  and  describe  it  here  briefly.  By  picking  a 
box  containing  the  unstable  saddle  with  part  of  its  stable  and 
unstable  manifolds,  we  can  determine  initial  conditions  that 
will  generate  trajectories  remaining  in  the  box  for  a  large 
number  of  iterations.  We  then  eliminate  any  points  converg¬ 
ing  to  an  attractor.  The  initial  conditions  remaining  in  the 
punctured  box  approximate  the  imion  of  the  stable  mani¬ 
folds,  while  the  last  point  of  the  trajectory  that  remains  in  the 
box  approximates  the  unstable  manifolds.  This  algorithm 
was  used  to  generate  the  stable  and  unstable  manifolds  in 
Fig.  8. 

As  shown  in  Fig.  8,  the  two-dimensional  stable  manifolds 
of  the  period-three  saddle  orbit  form  the  basin  bmmdary  be¬ 
tween  the  period-three  basin  and  period-four  basin.  The  one¬ 
dimensional  unstable  manifolds  approach  the  period-four  or¬ 
bit,  intersecting  the  two-dimensional  stable  manifolds  of  the 
period-two  and  period-one  saddles.  This  forms  a  forward 
connection  of  a  heteroclinic  tangle  in  R^.  There  are  no  re¬ 
verse  connections  or  intersections  of  the  stable  manifolds  of 
the  period-three  saddle  orbit,  which  would  be  necessary  for 
fully  developed  chaos. 

By  adding  stochastic  perturbations  with  a  large  enough 
standard  deviation,  random  trajectories  frequently  escape 
their  asymptotic  limit  toward  one  of  the  attracting  periodic 
orbits  and  visit  the  other.  In  contrast  to  basin  hopping,  the 
trajectories  spenda  significant  amount  of  time  in  between  the 
two  attractors,  near  the  forward  coimections  of  the  hetero- 
clinic  tangle.  Essentially,  short  visits  to  the  other  basin  act 
like  a  reverse  connection,  completing  the  tangle  and  enabling 
chaoslike  behavior.  Therefore,  the  trajectory  follows  the  cha¬ 
oslike  dynamics  in  the  time  spent  in  between  the  two  attrac¬ 
tors.  As  the  standard  deviation  of  the  noise  is  increased,  this 
reverse  jump  occurs  more  frequently  and  more  time  is  spent 
in  between  the  attractors.  These  events  can  be  identified  by 
bursting,  and  the  chaoticlike  behavior  is  captured  by  both  the 
Lyapunov  exponent  and  burst-rate  statistics. 

What  we  wish  to  identify  here  is  where  the  noise  facili¬ 
tates  the  reverse  jump  in  phase  space  and  provide  evidence 
that  the  phenomenon  is  similar  to  a  heteroclinic  tangle.  We 
begin  by  analyzing  the  time-series  data.  In  Fig.  8,  the  point 
before-the  trajectory  switches  basins  (defined  in  the  noise- 
free  case)  is  recorded.  It  is  clear  that  the  jumps  occur  fre¬ 
quently  in  three  regions  near  the  imstable  period-three 
saddles.  To  quantify  these  regions,  we  calculate  the  Galerkin 
transport  matrix. 

The  Galerkin  transport  matrix  can  be  used  as  a  tool  to 
identify  transport  between  the  original  basins  as  a  function  of 
the  standard  deviation  of  the  noise  added  to  the  system  [18]. 
Let  v(x)  be  the  distribution  of  the  random  variable  17, 

Kx)  =  e-(*''^'‘*^'2/^(2rr)3det(S) .  (6) 

As  a  spatial  approximation,  we  use  the  SFPO  in  the  form 
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FIG.  8.  (Color  online)  The  period-three  basin 
is  denoted  by  stripes  and  the  rest  of  the  space  is 
the  period-four  basin.  The  union  of  the  stable 
manifolds  in  the  phase  space  is  approximated  by 
small  dots  in  the  period-four  basin.  This  includes 
the  stable  manifolds  for  the  period-one  and 
period-two  orbits.  The  union  of  the  unstable 
manifolds  is  the  solid  curve  with  the  period-three 
stars  at  the  ends.  The  boundary  between  the 
period-three  and  period-four  basins  is  formed  by 
the  stable  manifold  of  the  period-three  saddle 
(squares).  Notice  how  the  unstable  manifolds  of 
the  period-three  saddle  intersect  the  stable  mani¬ 
folds  inside  the  period-four  basin,  forming  the 
forward  connections  of  the  heteroclinic  tangle. 
Superimposed  is  the  time-series  approximation  of 
the  flux  from  the  period-four  to  the  period-three 
basin  in  large  dots. 


PfJ.p(i0]=  Kix,y)piy)dy,  (7) 

where  the  stochastic  kernel  describing  the  PDF  of  the  noise 
perturbation  is  K{x,y)  =  i>(x-F{y)).  Assuming  a  nonzero  ex¬ 
ternal  noise  is  added  in  each  component,  Eq.  (7)  becomes 


Note  that  although  it  is  possible  to  let  any  of  the  standard 
deviations  tend  to  zero  in  the  SFPO  where  the  kernel  limits 
to  a  delta  function,  it  is  more  realistic  to  approximate  the 
zeros  by  very  small  values.  This  is  due  to  the  fact  that  the 
experiment  is  always  perturbed  by  small  noise.  Since  we 
require  a  finite  dimension  for  computation,  we  cover  the 
phase  space  with  N  disjoint  boxes  5,-  and  choose  a  set  of 
basis  functions  to  be  the  family  of  characteristic  functions 


(piix)  = 


f  1  if  X  e  Bj 
[o  if  X  $  Bj. 


(9) 


In  principle,  any  set  of  basis  fimctions  of  can  be  used,  but 
we  use  characteristic  functions  to  help  us  locate  spatial  trans¬ 
port,  as  was  motivated  historically  by  Ulam’s  method.  The 
approximation  of  the  Frobenius-Perron  operator  projects  to  a 
VX  A  matrix,  called  the  Galerkin  transport  matrix  (GTM), 


PpJ_4>,{x)'](f>j{x)dx, 


(10) 


The  GTM  describes  the  mass  flow  from  one  box  to  another 
over  one  iteration.  That  is,  the  entry  for  Ajj  approximates  the 
percentage  of  box  i  that  iterates  to  box  j  under  the  stochastic 
map.  Then  partition  the  boxes  according  to  their  basin  and 
reorder  the  GTM  by  similarity  transformations  to  reflect  that 


partition.  In  theory,  the  case  with  no  noise  will  result  in  a 
block-diagonal  matrix,  reflecting  dynamics  in  the  disjoint  ba¬ 
sins.  Under  stochastic  perturbations,  the  GTM  approximates 
three  things:  (i)  the  off-diagonal  blocks  indicate  where  the 
transport  between  basins  occurs — this  is  the  mass  flux  (or 
simply  flux),  (ii)  the  dominant  eigenvector  having  eigenvalue 
unity  approximates  the  PDF,  and  (iii)  by  weighting  the  mass 
flux  by  the  PDF,  we  pinpoint  regions  in  phase  space  that 
have  the  greatest  probability  of  leakage  into  another  basin — 
this  is  the  area  flux.  See  [19,18]  for  details. 

We  show  the  GTM  approximation  of  the  PDF  in  Fig.  9. 
Since  the  noise  distribution  is  assumed  to  be  normal,  it  is 
expected  that  the  PDF  has  nonzero  entries  everywhere  in 
phase  space.  However,  many  of  these  values  are  sufficiently 
small  so  that  when  added  to  unity,  they  make  zero  contribu¬ 
tion  due  to  the  fact  they  are  below  machine  error.  Therefore, 


FIG.  9.  An  approximation  by  the  GTM  of  the  PDF  when  o-j 
=0.16.  The  squares  represent  the  stable  period-four  orbit.  The 
circles  represent  the  stable  period-three  orbit.  The  darker  shades 
indicate  regions  with  the  highest  probability.  Notice  that  they  occur 
near  the  stable  periodic  orbits,  but  there  is  structure  connecting 
these  regions  called  the  stochastic  chaotic  saddle.  As  cr  increases, 
orbits  spend  more  time  on  the  chaotic  saddle,  indicating  the  in¬ 
creased  frequency  of  bursting  dynamics. 
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FIG.  10.  (Color  online)  Approximation  of  the  transport  from  the 
GTM  when  cri=0.16.  (a)  shows  regions  of  most  active  transport 
from  the  period-three  to  the  period-four  basin  in  large  dots,  (b) 
shows  regions  of  most  active  transport  from  the  period-four  to  the 
period-three  basin  in  large  dots.  Notice  that  these  regions  occur  near 
the  period-three  saddle  orbit  represented  by  the  squares.  In  (a),  the 
union  the  stable  manifolds  is  displayed  in  layered  sheets  and  un¬ 
stable  manifolds  form  the  one-dimensional  curve  with  the  period- 
three  points  (stars)  at  the  ends.  They  are  approximated  by  the  box 
algorithm  from  [13].  The  stable  manifolds  of  the  period-three 
saddle  form  the  basin  boundary  between  the  two  basins.  In  (b),  the 
basin  of  the  period-three  orbit  is  represented  by  the  small  dots.  Also 
shown  in  both  (a)  and  (b)  are  the  stable  period-four  orbit  (stars),  the 
unstable  period-two  orbit  (triangles),  and  the  unstable  period-one 
(circle). 

we  choose  a  numerical  threshold  of  machine  precision  as  a 
lower  botmd  and  replace  all  smaller  values  to  zero  in  the 
PDF.  Notice  that  as  we  add  stochastic  perturbations  to  the 
system,  the  most  frequently  visited  regions  lie  near  the  un¬ 
stable  manifolds  of  the  period-two  and  period-one  saddle 
orbits  from  the  noiseless  case. 

We  observe  that  as  the  standard  deviation  is  increased,  the 
PDF  spreads  and  crosses  into  the  period-three  basin,  and 
through  the  stable  manifold  of  the  period-three  saddle.  This 
is  evidence  that  there  is  the  topology  for  a  trajectory  to  emu¬ 
late  chaoslike  behavior.  We  now  will  verify  the  fact  that  the 
trajectory  actually  uses  these  regions  for  transport.  This  is 
supported  by  the  area  flux,  which  is  shown  in  Fig.  10.  The 


FIG.  II.  Bursting  time  series  from  the  experiment  (a)  and  model 
(b).  The  dashed  line  indicates  the  threshold  used  to  determine  a 
burst. 


regions  where  the  trajectory  is  most  likely  to  switch  basins 
are  found  by  multiplying  the  mass  flux  by  the  associated 
PDF  value  for  that  region  of  phase  space.  Notice  the  agree¬ 
ment  between  the  transport  region  predicted  by  the  time  se¬ 
ries  in  Fig.  8  and  the  area  flux  from  the  period-four  to  the 
period-three  basin. 


VI.  DISCUSSION 

Dynamics  with  noise  is  always  present  in  experiments  at 
least  at  some  level.  In  many  cases,  noise  is  sufficiently  small 
so  that  its  role  is  ignorable  with  respect  to  the  underlying 
determinism.  However,  even  relatively  low-amplitude  noise 
may  play  a  significant  role  in  which  the  dynamics  takes  on  a 
qualitative  change  that  is  different  from  the  deterministic 
structure.  In  the  physical  example  presented  here,  we  have 
examined  an  experiment  where  noise  has  been  injected  into  a 
modulated  laser.  The  amplitude  of  the  noise  was  adjusted, 
and  the  laser  was  seen  to  go  from  stochastically  perturbed 
periodic  behavior  to  one  of  stochastic-induced  chaoslike  dy¬ 
namics.  Because  discrete  control  of  the  random  noise  ampli¬ 
tude  could  be  achieved,  the  system  was  therefore  analyzable 
by  a  discrete-map  approach,  thereby  revealing  explicitly  the 
interaction  of  noise  and  the  underlying  deterministic  topol¬ 
ogy- 

In  conjunction  with  the  laser  experiment,  we  have  exam¬ 
ined  a  quantitative  model  with  additive  noise  in  the  intensity. 
Both  exhibit  similar  bursting  behavior,  as  shown  in  the  time 
series  data  in  Fig.  11.  Although  the  topology  of  the  experi¬ 
mental  dynamics  is  difficult  to  ascertain,  the  quantitative  na¬ 
ture  of  the  model  does  allow  an  in-depth  view  of  the  imder- 
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FIG.  12.  (Color  online)  The  return  map  of  the  fixed  points  cor¬ 
responding  to  the  periodic  orbits  and  their  projected  manifolds. 
(The  X  axis  has  been  plotted  logarithmically  to  show  more  detail.) 
The  unstable  manifolds  form  the  dark  solid  curve  in  the  middle, 
while  the  stable  manifolds  are  approximated  by  smaller  points.  In 
addition,  the  projected  regions  of  transport  are  overlaid  in  large 
dots.  Notice  that  the  transport  between  the  two  basins  predicted  by 
the  GTM  in  Sec.  V  lies  close  to  the  period-three  saddles  (triangles) 
on  the  stable  manifold.  The  other  periodic  orbits  are  labeled  as 
follows;  period  two,  stars;  period  four,  squares;  period  one,  large 
dot;  period-three  node,  circles;  period-three  saddle,  triangle.  The 
value  of  the  standard  deviation  used  was  0-1=0.04. 

lying  topology  and  its  relation  to  noise.  In  the  absence  of 
noise,  the  topology  of  the  system  was  determined,  and  the 
structure  of  the  stable  and  unstable  manifolds  was  computed 
in  a  prechaotic  regime.  When  noise  is  added,  the  stmcture  of 
the  topology  interacts  with  the  stochastic  fluctuations  in  such 
a  way  to  induce  chaoslike  behavior,  which  is  the  emergent 
structure  observed  in  both  theory  and  experiment  as  shown 
in  Fig.  4.  The  stochastic  dynamics  is  the  union  of  local  sto¬ 
chastic  dynamics  within  each  basin  and  the  dynamics  near  a 
partially  formed  chaotic  saddle.  For  sufficiently  large  noise 
amplitudes,  local  instability  near  the  manifold  structure  con¬ 
tributes  to  the  time-  and  space-averaged  linear  variation  so 
that  the  Lyapunov  exponent  becomes  positive,  which  we  take 
as  criteria  for  stochastic  bifurcation  as  defined  in  [15],  and 
exemplified  in  [2]. 

In  tying  together  the  dynamics  from  the  model  and  ex¬ 
periment,  we  can  project  the  phase  portrait  of  the  transport 
and  manifold  structure  to  a  lower  dimensional  return  map,  as 
we  did  in  Sec.  IV.  In  Fig.  12,  we  have  depicted  the  fixed 
points  that  correspond  to  the  unstable  periodic  orbits  and 
their  manifolds.  Notice  that  because  we  used  the  box  algo¬ 
rithm  of  [13],  the  manifolds  are  not  grown  from  a  saddle,  but 
reflect  the  union  of  all  such  manifolds  in  the  region  we  con¬ 
sidered.  The  stable  manifold  (in  black)  corresponds  to  the 
basin  boundary  in  the  original  phase  space  separating  the 
bistable  attractors  in  the  deterministic  case.  In  the  projection 
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in  Fig.  12,  we  can  see  where  the  maximum  probability  of 
transport  from  one  basin  to  another  occurs.  Notice  that  much 
of  it  takes  place  near  the  basin  boundary  saddle  points.  That 
is,  the  period-three  saddle  stable  manifold,  which  forms  the 
basin  boundary,  intersects  the  regions  of  maximal  probability 
transport.  The  stochastic  dynamics  fluctuates  until  it  comes 
near  the  basin  boundary,  at  which  point  it  is  attracted  to  the 
saddle  point.  Noise  then  takes  on  a  dominant  role,  where  the 
intensity  is  either  pushed  across  the  basin  boundary,  or  re¬ 
mains  in  the  same  basin.  The  unstable  manifold  then  domi¬ 
nates  the  noise,  by  pushing  the  dynamics  further  into  the 
respective  basin. 

Our  model  of  the  class-5  laser  includes  more  physics  in 
the  problem,  which  in  turn,  leads  to  a  more  interesting  class 
of  dynamical  behavior.  Most  laser  models  consist  of  just  two 
population  levels,  describing  the  change  in  the  population 
inversion  and  intensity.  The  advantage  of  such  a  model  is 
that  it  may  be  studied  in  the  plane,  having  one-dimensional 
stable  and  unstable  manifolds.  The  resulting  stochastic  analy¬ 
sis  could  be  compared  directly  to  manifolds  that  were  built 
from  curve-following  methods  in  the  plane.  In  contrast,  the 
model  considered  here  is  based  on  a  four-level  model,  which 
agrees  quantitatively  with  the  experiment  over  a  large  range 
of  values.  The  model  requires  two  main  levels  and  two  rota¬ 
tional  levels,  resulting  in  a  five-dimensional  system  of  differ¬ 
ential  equations.  Approximating  the  relaxation  rates  of  the 
vibrational  states  by  their  average  allows  one  to  reduce  the 
model  to  the  current  three-dimensional  driven  case  [24].  The 
main  difference  here  is  that  the  stable  manifolds  are  no 
longer  one-dimensional.  (The  unstable  manifolds  are  one¬ 
dimensional,  however.)  Here,  two-dimensional  stable  mani¬ 
folds  are  pierced  by  one-dimensional  unstable  manifolds. 
Therefore,  regions  of  transient  behavior  may  wander  over  a 
greater  region  of  phase  space  in  both  the  deterministic  and 
stochastic  models,  offering  a  richer  set  of  dynamical  behav¬ 
ior  than  the  two-level  laser  model. 

One  of  the  main  conclusions  of  the  cizrrent  stochastic 
analysis  is  that  maximal  transport  from  one  basin  to  another 
may  not  occur  near  the  basin  boundary  saddles.  Similar  ex¬ 
amples  based  on  asymptotic  properties  of  problems  of  escape 
where  the  phenomenon  of  saddle  avoidance  occurs  can  be 
found  in  [28,29].  We  note  that  the  methods  used  here  not 
only  agree  with  the  previous  local  theories,  but  is  an  alterna¬ 
tive  to  describe  the  global  structure  of  the  transport  as  well 
[18]. 

In  general,  computing  stable  and  unstable  manifolds  is  a 
difficult  task,  compounded  here  by  the  fact  that  the  mani¬ 
folds  are  of  different  dimensions.  The  technique  used  in  this 
paper  cannot  grow  the  manifolds  from  a  given  saddle.  There¬ 
fore,  the  global  analysis  of  the  four-level  laser  system  lacks 
some  of  the  precision  of  the  two-level  system.  On  the  other 
hand,  the  SFPO  tool  does  not  require  the  manifold  construc¬ 
tions  a  priori.  Rather,  the  transport  requires  a  partition  of  the 
basins  of  attraction  in  the  zero  noise  case.  When  noise  is 
added,  the  phase  space  is  reconstructed  in  terms  of  transport 
across  the  referenced  basin  boundaries,  and  thus  must  con¬ 
tain  components  of  the  basin  manifolds,  regardless  of  their 
dimension. 

In  terms  of  model  development,  much  work  on  noise  had 
been  done  on  maps  in  the  plane,  or  two-dimensional  flows. 
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In  this  paper,  we  have  presented  a  stochastic  experiment  as 
well  as  a  quantitative  model  that  simulates  the  stochastic 
dynamics.  The  model  is  itself  a  reduction  of  a  previously 
more  complicated  model  [24],  but  nonetheless,  captures  the 
relevant  features  of  the  stochastic  dynamics.  In  particular,  it 
captures  the  interaction  of  the  stochastic  dynamics  with  the 
underlying  topology  of  the  model.  This  quantifies  the  smooth 
change  of  regular  stochastic  behavior  to  a  bursting  type  of 
behavior  between  basins,  which  appears  to  be  chaoslike  in 


that  there  exists  a  smooth  transition  from  negative  to  positive 
Lyapunov  exponents. 
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